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c A model of a ringdown signal with noise is created to mimic 

c experimental data obtained from 

c a photodector and preamplifier applied to an A to D convertor. 

c It is then fit with the new algorithm that subtracts a line 

c that is fit to the background before applying the standard 

c algorithm and compares that fit with one that simply uses a 

c constant for baseline subtraction. We use exponential decays 

c of tau=80 microseconds, choose a peak signal, add gaussian noise, 

c add a dc offset then add sixty Hz noise. As in real experiments 

c the data is sampled every microsecond. 

c Then use one iteration of nonlinear fit and compare. 

c ( Note: Fortran 77, with a few exceptions, is used. One exception 

c which is machine specific is the open statement. As an example open 

c (39,file=' ' , status= ' new ' ) . The file name which is left 

c blank allows its insertion at run time.) The text that follows an "!" 

c is a comment and not part of the source code. 

c 

program model 
implicit none 

integer iy , i , max , k , j , j max , j j , max2 , flag 

parameter (k=2000,jj =500) ! k= number of data pts sampled in 
c microseconds, jj= number of repeats to form ensemble 

real*8 s(k) ,av(jj) ,amp,offst,ac,avtau0,x,ava,d,ss0,ssl ,y 
real*8 s2(k) ,s3(k) ,m,b,chisq(jj) , sural , sum2 , sum3 , sumt , sumtO , bb , mmm 
real*8 sum4 , sum5 , t , a (j j ) , sigma , tauO (j j ) , sigmatO , dum , z , sigmaa , gamp 
real* 4 urand,ss 
external function urand 

write (*,*) 'pick an integer for random number generator seed' 
read(*,*) iy 

flag=0 ! flag = 0 uses new algorithm, flag = 1 uses old 

c get the data choose added 60hz noise etc. 

print* ,' Enter the Gaussian white noise sdev in milivolts' 
read(*,*) gamp 

print* ,' Enter signal amplitude' 
read(*,*) amp 

print*, 'Enter 60 Hz rms level' 
read(*,*) ac 

print* , ' Enter the DC offset voltage ' 
read (* , *) of fst 

t=80.d0 ! decay time 

write (*,*) ' Name a result file' 
open (39 , f ile= ' ' , status= ' new ' ) 
write (39,*) The chosen values are:' 

write (39.*) ' random seed =' ,iy, ' , amp =' .real (amp) , '60Hz-rms =' , 
*real (ac) , ' offset = ' , real (of fst) , ' tau = ' , real (t) 

write (39,*)' tau Amp num points sigma reduced chisq 

*sigma-tau sigma-amp element num' 

write (39,*)' ' 

ac=dsqrt(2.d0)*ac ! convert to amplitude, we entered rms value 
550 continue 

500 continue 



j = l+j 
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call getdata (offst,t, amp, ac, gamp, iy,s(l)) [produce the model data 
want to fit the record s to the exponential, first using a weighted log transformatl 



c get average value for the last 100 pts of each decay for base line determination, 
c For these ponts signal has fallen many orders of magnitude below the noise. 
av(j)=0.d0 
do 35 i=k-99,k 
av(j)=av(j)+s(i) 
35 continue 

av(j)=av(j)/100.d0 



c subtract base line from the data and find sigma from last 100 pts 

sigma=0 . dO 

do 45 i=l,k 
s2(i)=s(i)-av(j) 

if (i.ge.k-99)then 

s igma=s igma+ (s ( i) - av (j ) ) * * 2 

endif 

45 continue 

sigma=sigma/(100.d0-l.d0) ! calculate sigma from baseline fluctuation 
sigma=sqrt (sigma) 



c check for values of data below 1 sigma to obtain max , highest numbered positve 
c non zero point 



i=l 

while (s2 (i) . ge . 1 . d0*sigma) 

max=i- 1 

i=i+l 

enddo 



if (flag.eq.l) goto 560 ! skip line fitting and subtraction if done once 



c fit straight line to data and subtract line from function rather than just a 

c constant 

if (k-2*max.le.200) then 
max2=200 

else ! begin data for line at twice the time that exponential reaches noise 
max2=2*max ! at this point signal is as far below noise as it was above 
!the noise initially, 
endif ! but use no less than 200 pts for fit. 



algorithm for straight line fit x -> i, y -> s(i) y=(mmm)x+bb 

suml=0.d0 
sum2=0 . dO 
sum3=0 . dO 
sum4=0 . dO 
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do 52 i=max2,k 

suml=suml+s(i) ! sum(y) 

sum2=sum2+i ! sum(x) 

sum3=sum3+i**2 ! sum (x A 2) 

sum4=sum4+i*s(i) ! sum (xy) 

52 continue 
c Npt's=(k- max2+l) 

d=(k- max2+l) *sum3-sum2**2 

if (d.eq.O.dO) then 

print* ,' singular solution' ! make sure program doesn't crash if no solution 

pause 

end if 

mmm=((k- max2+l) *sum4-suml*sum2) /d ! slope 

bb=(sum3*suml-sum2*sum4)/d ! intercept 

c 

c now subtract line from the initial data 

do 54 i=l,k 

s2(i)=s(i)-bb-i*mmm ! line y = mmmx +bb 
54 continue 



while (s2(i) .ge. 1 .d0*sigma) ! find maximum value of i for log fit 

max=i- 1 
i=i+l 
enddo 

c 

560 continue 

c now do log transform and calculate slope, and intercept of line with weighted 

c average 

c 

c 

do 60 i=l,max 

s3(i)=dlog(s2(i)) ! convert to log of function 
60 continue 



sum 1=0. dO 
sum2=0 . dO 
sum3=0 . dO 
sum4=0 . dO 
sum5=0 . dO 

= prelog y A 2/sigma A 2 = s2(i) A 2/sigma A 2 
common factor sigma A 2 taken out from all terms 
! sum (w*x) 
! sum (w*y) 
! sum (x A 2*w) 
! sum w*x*y 



do 65 i=l,max ! w=weight 

suml=suml+s2(i)**2 ! sum( w) 
sum2=sum2+i*s2 (i) **2 
sum3=sum3+s2 (i) * *2*s3 (i) 
sum4=sum4+i**2*s2(i) **2 
sum5=sum5+i*s2 (i) **2*s3 (i) 
65 continue 



m= (sum 1 * sum5 - sum3 *sum2) / (sum 1 * sum4 - sum2 * * 2) 
b=(sum3/suml) - (m*sum2/suml) 



c Also calculate reduced average chisq and sigma tau, see Bevington 

chisq (j)=0.d0 
do 75 i=l,max 

chisq (j ) = chisq (j ) + (s2 (i) -dexp (b) *dexp (m*i) ) * *2/ ( (max- 2) * (sigma** 2) ) 
75 continue 



tau0(j)=-l .d0/m ! convert from rate to tau 
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a (j ) =dexp (b) 

x=sigma*t**2*dsqrt(suml/(suml*sum4-sum2**2)) ! from Bevington ,sigma t 
y =sigma*a(j)*dsqrt(sum4/(suml*sum4-sum2**2)) ! page 108-9, sigma amp 

write (39,*) real (tauO (j ) ) , real (a (J ) ) , ' ',max,' 
*', real (sigma) , ' ' , real (chisq (j ) ) , ' ' ,real(x) , ' \real(y),' ' , j 



jmax=j 

if (j.ge.jj) goto 100 

goto 500 ! get new data to develop new ringdown for the ensembie 

100 continue 



c calculate ensemble averages 

sumt0=0 . dO 
dum=0 . dO 
sumt=0 . dO 

do 120 j=l , jmax !find average tau's, av's , a's, and chisq's 
sumt=sumt+a(j) 
sumt0=sumt0+tau0 (j ) 
dum=dum+chisq (j ) 
120 continue 

dum=dum/jmax ! average chisq 

avtau0=sumt0/j max 
ava=sumt/j max 

sigmaa=0.d0 
sigmat0=0.d0 

do 130 j = l, jmax !find sigma tau and sigma amp 

sigmat0=sigmat0+ (tau0(j ) -avtauO) **2 
sigmaa=sigmaa+ (ava-a (j ) ) * *2 
130 continue 

s igmaa=dsqrt (s igmaa/ (j max - 1 ) ) 
sigmat0=dsqrt (sigmatO/ (j max- 1) ) 
if (flag.eq.0) ss0=sigmat0 
if (flag.eq. 1) ssl=sigmat0 
if (flag.eq. 1) ss=real(ssl/ss0) 

write (39,*) 'average tau = 1 , real (avtauO) , ' sigma tau = ', 
*real(sigmat0) , 'average red chisq = ' ,real(dum) 

write (39,*) 'average amp = ' , real (ava) , ' sigma amp = ' , real (s igmaa) 
if (flag.eq.0) then 

flag=l ! We now do it all over again using the simple constant baseline subtraction, 
write (39,*) ' ' 

write (39,*) 'These next results just subract a constant and not a line to fit 
write (39,*)' ' 

z=urand(iy+17) ! call a random number with a new seed to cause sequence to reset 
c when iy is used again 



c This insures that the data used the second time around will be identical to that used 

c the first fit. 

c In that way the comparison is exact between the two algorithms 



C 



goto 550 
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end if 

write (39,*) ' the squared ratio of sigma tau for the two data sets = ' , ss**2 
close (39) 

end 



subroutine getdata (of f st , t , amp , ac , gamp , iy , s) 

c 

c Given the parameters for the dc offset, decay time, signal amplitude, 60 hz. 

c amplitude, rms 

c random gaussian noise, and random number seed, returns an array S(k) that represents 

c the ringdown voltage. 

c 

c 

implicit none 
integer iy , istat , i , k 
parameter (k=2000) 
real*8 of fst , t , amp , ac , phi , gamp 
real *8 s (k) , g (k) , nO (k) , y (k) 
real* 4 gasdev, urand 
external function urand 
external function gasdev 

phi=dreal (urand (iy) *6 . 28) !the random phase of the 60Hz signal, a number between 0 and 2pi 
do 15 i=l,k 

g (i)=gamp*dreal (gasdev (iy) )! gaussian noise sd=l,avg value =0 ,for k micro sec 
y (i) =amp*dexp (-dreal (i) /t) + offst 

n0(i)=ac*dsin(0.000377d0*i+phi) ! sixty cycle with random phase for each ensemble member 
s(i)=g(i)+y(i)+nO(i) 
1 5 continue 
return 
end 



function gasdev (iy) 

c 

c generate normal gaussian random numbers sd=l avg =0 (From "Numerical Recipes" , Press et. al 

c Cambridge University Press 1986) 

c 

external function urand 

real* 4 urand , vl , v2 , gasdev , gset , fac , r 

integer *4 iset.iy 
data iset/0/ 
if ( iset . eq . 0) then 
1 vl=2 . *urand(iy) -1 . 

v2=2.*urand(iy)-l. 

r=vl**2+v2**2 

if (r.ge. 1 .)go to 1 

f ac=sqrt (- 2 . * log (r) /r) 

gset=vl*fac 

gasdev=v2*fac 

iset=l 
else 

gasdev=gset 

iset=0 
endif 
return 
end 
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real function urand(iy) 

c 

c Finds a psuedo random number. 

c the integer iy should be initialized to an arbitrary integer prior to the fist caii 
c to urand. the calling program should not alter the value of iy between subsequent 
c calls to urand. values are returned between 0 and 1. 



c 



integer ia, ic, itwo,m2 ,m,mic 
double precision halfm 
real s 

double precision datan.dsqrt 
data m2/0/ , itwo/2/ 
if(m2.ne.O) go to 20 



c if first reentry compute machine integer word length 



c 

m=l 

10 m2=m 



m=itwo*m2 

if (m.gt.m2) go to 10 
halfm=m2 



ia=8* idint (halfm*datan ( 1 . dO) /8 . dO) +5 

ic=2* idint (halfm* (0 . 5d0 -dsqrt (3 . dO) /6 . dO) ) + 1 

mic=(m2-ic) +m2 

c 

c s is a scale factor for converting to floating point 
c 

s=0 . 5/halfm 

c 

c compute next random number 
c 

20 iy=iy*ia 

if (iy.gt.mic) iy=(iy-m2) -m2 
iy=iy+ic 

if (iy/2.gt.m2) iy=(iy-m2) -m2 

if (iy . It . 0) iy= (iy+m2) +m2 

urand=f loat (iy) *s 

return 

end 
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This is the compiled results using 500 independent ringdowns with the rms 60 Hz = to 

the eaussian rms white noise. As you can see, as indicated at the bottom of the table, the 

squared ratio of the ensemble sigma tau is 19.77. That ratio shows that to achieve a given 

precision in tau, one needs 19.77 times as much averaging using the standard algorithm 

comparted with the new one. We have repeated this run 10 times to develop a statistic on the 

improvement for this set of parameters and the result was that the average improvement was 

18 98 with a sdev = 1.42. As you might expect, if the 60 Hz signal gets smaller, the 

improvement also gets smaller. However, the signifcance of the result is that the time taken 

to make a measurement of a given precision is reduced by that ratio. 

(This is the output of the source code with one item that I added subsequently , 

the "gauss rms = 3.5). 



The chosen values are: 
random seed =1 ,amp = 800.000 60Hz-rms 



3.50000 offset 



tau 




Amp 


num points 


sigma 


reduced chisq 




JO 1 c 


799. 


880 


377 


3. 


41363 


1 . 


06457 


oU . 




800. 


545 


370 


3. 


47633 


1 . 


11195 


ou . 




797. 


733 


357 


3. 


56602 


0. 


901096 


70, 




800. 


223 


364 


3. 


47949 


0. 


935880 


on 
oU . 


97^7 


799. 


792 


343 


3. 


85445 


0. 


910839 


7Q 


71 


800. 


013 


356 


3. 


64096 


0. 


903653 


on 
oU . 


1 H7Q 


799. 


414 


328 


3. 


66295 


0. 


902719 


on 
oU . 


4UoD 


798. 


814 


369 


3. 


40928 


1 . 


09774 




yo^u 


798. 


997 


354 


3. 


46550 


0. 


919112 


vn 

7y . 


94DU 


799. 


908 


339 


3 


58053 


0. 


859031 


o f~\ 
OU . 


5 /DO 


799. 


571 


369 


3 


49413 


0. 


977574 


on 




799. 


982 


367 


3 . 


42647 


1 . 


11993 


80 . 


lUbo 


799. 


776 


364 


3 


68315 


0. 


979039 


7Q 

/ y . 


S2S29R 


797. 


612 


355 


3. 


67865 


0. 


896638 


on 


iy4U 


800. 


994 


342 


3 


32195 


1 , 


16005 


on 


1 1 QO 

i lay 


799. 


735 


368 


3. 


6681 1 


0 , 


845112 


oU . 


99A.9. 


797. 


985 


370 


3. 


66111 


0. 


990751 


80. 


2534 


798. 


529 


337 


3. 


09522 


1. 


33661 


80. 


6539 


798. 


334 


377 


3. 


34951 


1 . 


,02779 


80. 


0661 


798. 


940 


352 


3. 


36445 


1 


,05537 


80, 


1926 


799. 


261 


344 


3. 


70351 


0 


.957610 


79 


.7549 


800, 


.642 


390 


3 


.43053 


0 


.941505 


79. 


.9710 


799, 


.103 


360 


3 


. 16224 


1 


. 03390 


80 


.2389 


799 


. 133 


352 


3 


. 70208 


0 


.869751 


80 


.3320 


799 


.336 


362 


3 


.48165 


0 


.898230 


80 


.0478 


799 


.761 


357 


3 


.55346 


0 


.874637 


80 


.1229 


801 


.448 


358 


3 


. 57680 


0 


.925075 


79 


.8343 


799 


.408 


362 


3 


.52717 


0 


.915940 


80 


.4072 


800 


.646 


345 


3 


.23790 


1 


.24147 


80 


.2135 


798 


.802 


362 


3 


.46311 


1 


.02765 


79 


.9797 


798 


.023 


347 


3 


.59251 


0 


.993194 


80 


.0015 


800 


.016 


335 


3 


. 17675 


1 


.31002 


80 


.2739 


799 


.292 


370 


2 


. 80468 


1 


.36614 


80 


.2325 


799 


.895 


355 


3 


. 18023 


1 


.28376 


80 


.1239 


799 


.947 


357 


4 


.01733 


0 


. 783859 


80 


.3681 


798 


.786 


359 


3 


.29321 


1 


.11531 


80 


.1854 


797 


.797 


345 


3 


.85224 


0 


.710844 


80 


.4123 


799 


.192 


352 


3 


.47029 


1 


.01612 


80 


.0014 


799 


.991 


384 


3 


. 18610 


1 


.05959 


80 


.1084 


799 


.812 


348 


3 


.39072 


1 


.05840 


80 


. 1952 


800 


.702 


354 


3 


.52682 


1 


.07649 


80 


.3925 


798 


.095 


353 


3 


.50967 


0 


. 866208 


80 


.0717 


798 


.897 


354 


3 


.22297 


1 


.40809 



10.0000 tau = 
sigma-tau 



80.0000 gauss rms =3.5 
sigma-amp element num 



0. 108180 
0.110807 

0.113038 

0.112026 

0. 122796 

0. 117230 

0.117561 
0.108281 

0.111112 

0.115175 

0.110215 
0. 108919 

0. 117587 

0. 118441 
0. 105765 

0. 117201 

0. 116966 
9.911115E-02 
0.105757 
0. 107878 

0. 118209 

0. 109912 
0. 101494 

0. 118023 

0. 110457 

0. 113643 

0. 113715 

0. 113437 
0.102770 
0.110439 

0. 115368 
0.101899 
8.940282E-02 
0.101511 

0. 128420 
0. 104797 

0. 123417 
0. 110295 
0.101856 
0.108519 
0.112333 

0. 111681 
0. 103220 



0.771510 
0.787399 
0 . 805365 
0.791412 
0.872772 
0.827762 
0.831779 
0.770851 
0.786289 
0.814045 
0 . 788068 
0.775385 
0.834578 
0.835436 
0.752381 
0.831674 
0.829375 
0.702039 
0.755472 
0.763572 
0.838947 
0.778532 
0.717973 
0.838264 
0.786854 
0.805859 
0.809729 
0.801489 
0.732525 
0 . 784080 

0.815205 
0.721247 

0.635305 
0.720895 

0.911108 
0.745206 

0.873547 
0.784959 
0.722510 
0.769263 
0.798850 

0.793859 
0.730957 



1 

2 

3 
4 
5 
6 
7 

8 

9 

10 

11 

12 
13 
14 

15 
16 
17 

18 

19 

20 
21 
22 

23 
24 
25 
26 
27 
28 

29 

30 
31 

32 
33 

34 
35 

36 
37 

38 

39 

40 

41 
42 

43 
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80 


.0652 


799 


.503 


361 


3 


01224 


1 


.36238 


9.651958E-02 


0.683667 


44 


79 


.9924 


799 


581 


348 


3 


52276 


0 


.902304 


0. 113262 




0.800741 


45 


80 


.2594 


799 


.905 


360 


3 


.80231 


0 


.832782 


0. 120938 




0.860455 


46 


79 


.8885 


798 


848 


360 


3 


34321 


1 


.14142 


0.107473 




0.759256 


47 


79 


7865 


799 


934 


362 


3 


.80141 


0 


.846560 


0.122173 




0.863599 


48 


80 


.0431 


798 


649 


356 


3 


69423 


0 


.826237 


0. 118471 




0.838282 


49 


80 


.1002 


798 


542 


357 


3 


65684 


1 


.00040 


0. 117087 




0.829213 


50 


80 


.1506 


799 


860 


350 


3 


42042 


1 


.25200 


0.109397 




0.775831 


51 


79 


8712 


799 


471 


368 


3 


43632 


0 


.969244 


0. 109895 




0 . 778789 


52 


80 


2386 


800 


389 


373 


3 


52289 


1 


.00298 


0. 111996 




0.797288 


53 


80 


.6664 


800 


708 


376 


3 


45900 


0 


.982347 


0. 108635 




0.779197 


54 


80 


2141 


799 


533 


354 


3 


50310 


1 


00037 


0. 111651 




0.793270 


55 


80 


1684 


799 


886 


376 


3 


57510 


0 


.800629 


0. 113680 




0.808862 


56 


79 


8979 


798 


265 


379 


3 


53829 


0 


979120 


0.113534 




0.802698 


57 


79 


8790 


801 


421 


328 


3 


56524 


1 


12472 


0. 114764 




0.811164 


58 


80 


0871 


800 


029 


354 


3 


43421 


0 


927179 


0 . 109803 




0.778901 


59 


80 


2119 


798 


231 


367 


3 


36987 


1 


23473 


0. 107429 




0.762613 


60 


79 


8350 


798 


362 


366 


3 


41828 


0 


922450 


0. 109908 




0.776138 


61 


79 


7611 


798 


936 


358 


3 


38721 


1 


12772 


0. 109312 




0.770527 


62 


80 


0349 


800 


326 


341 


3 


67013 


0 


884792 


0.117623 




0.833229 


63 


80 


2091 


798 


672 


350 


3 


45391 


1 


06840 


0. 110298 




0 . 782376 


64 


80 


1261 


798 


263 


354 


3 


90308 


0 


721050 


0. 124728 




0.884100 


65 


79 


8974 


797 


609 


363 


3 


26460 


1 


20106 


0.105494 




0 . 742836 


66 


80 


2315 


800 


633 


383 


3 


38483 


0 


997836 


0. 107483 




0.765750 


67 


80 


2131 


800 


280 


345 


3 


77116 


0 


829861 


0. 120007 




0.853593 


68 


80 


3054 


800 


356 


357 


3 
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average tau = 80.1192 sigma tau = 0.197362 average red chisq = 1.01401 



average amp = 799.481 sigma amp = 0.814629 
These next results just subract a constant and not a line to fit 
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KJ . 1. XKJ\J tJ 1 


0.804778 


473 


i . iyj 10 


f) 1 1 3879 


0.793143 


474 


1 42168 


0 . 103587 


0.725598 


475 


1 D7S10 
1 . u / o ± u 


0 . 1 15057 


0 . 808890 


476 


1 20489 


0. 107867 


0.782929 


477 


1 .05364 


0. 111930 


0.780367 


478 


1.28707 


0.100545 


0.717853 


479 


0.988438 


0. 120240 


0.841224 


480 


1 .07175 


0. 114676 


0.802413 


481 


1 .29518 


0. 118836 


0.828453 


482 


0 . 824005 


0. 125253 


0.892442 


483 
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80 


6554 


800 


760 


349 


3 


94831 


0 


803779 


0.124404 


0.890676 


484 


79 


2522 


799 


895 


344 


3 


37891 


1 


16644 


0.110809 


0 . 774079 


485 


81 


4727 


799 


340 


386 


3 


21776 


1 


45687 


9.928988E-02 


0.719420 


486 


80 


3845 


798 


911 


343 


3 


26233 


1 


18222 


0.103644 


0.737549 


487 


79 


9768 


799 


515 


359 


3 


58749 


0 


924636 


0.115200 


0.814954 


488 


81 


2791 


799 


591 


397 


3 


66608 


1 


07796 


0.113404 


0.820533 


489 


81 


0803 


800 


646 


385 


3 


39049 


1 


15636 


0.105308 


0 . 760469 


490 


80 


0694 


800 


645 


348 


3 


26363 


1 


27675 


0.104382 


0.740426 


491 


81 


2979 


800 


465 


376 


3 


37366 


1 


12937 


0. 104534 


0 . 756077 


492 


80 


9854 


800 


178 


363 


3 


59515 


1 


12044 


0.112249 


0.807909 


493 


78 


9428 


799 


643 


355 


3 


61705 


0 


974909 


0. 118987 


0 . 829096 


494 


79 


0115 


799 


279 


333 


3 


81745 


0 


851071 


0. 125639 


0.874952 


495 


81 


4317 


798 


331 


406 


3 


66311 


0 


993203 


0.113023 


0.818486 


496 


81 


0048 


800 


026 


372 


3 


42858 


1 


18312 


0.106611 


0 . 768799 


497 


79 


0539 


798 


444 


325 


3 


43988 


1 


23858 


0.113469 


0.788877 


498 


79 


0262 


798 


720 


335 


3 


40804 


1 


12155 


0.112346 


0.781563 


499 


79 


7364 


797 


962 


339 


3 


60389 


1 


08535 


0.117128 


0.821964 


500 



average tau = 80.1321 sigma tau = 0.877687 average red chisq = 1.09949 

average amp = 799.267 sigma amp = 1.00279 

the squared ratio of sigma tau for the two data sets = 19.7766 



